%
% enhance(thresh)
%
% thresh = minimum value of reconstructed image that we won't zero out
%
% This function calculates the squared error between an original tomographic
% image (generated in MATLAB) and six reconstructed images.  The six
% squared error results are calculated using the following filters:
%
% r_sqerr = fft reconstruction using backproj.m
% p_sqerr = same as r_sqerr with negative values set to 0
% g_sqerr = same as p_sqerr with values lower than thresh set to 0
% m_sqerr = same as g_sqerr after being run through 3x3 median filter
% w_sqerr = same as g_sqerr after being run through 3x3 weiner filter
% h_sqerr = same as g_sqerr after being run Haar Wavlet edge sharpening
%           filter
%
% Developed by  Andy Doran
%               Rice University
% 		April, 1999
%

function enhance(thresh)

% Load the original image

pict = load('im1');
orig = pict.mat2;    % The string after the . denotes the
		     % image name.

% Load reconstructed image

rpict = load('rim1');
rimage = rpict.imag_250_100_250;  % The string after the . denotes the
				  % image name.

% Display squared error between original image and reconstructed image

error = orig - rimage;
r_sqerr = sum(sum(error.^2))

% Display squared error between original image and reconstructed image
% with all negative values zeroed.

pimage = zeroneg(rimage);
error = orig - pimage;
p_sqerr = sum(sum(error.^2))

% Display squared error between original image and reconstructed image
% with all negative values zeroed and all values less than thresh zeroed.

gimage = pimage;
inds = find(gimage < thresh);
gimage(inds) = 0;
error = orig - gimage;
g_sqerr = sum(sum(error.^2))

% Display squared error between original image and reconstructed image
% run through a 3x3 median filter

block = [3 3];
mimage = medfilt2(gimage,block);
error = orig - mimage;
m_sqerr = sum(sum(error.^2))

% Display squared error between original image and reconstructed image
% run through a 3x3 Weiner2 filter

block = [3 3];
wimage = wiener2(gimage,block);
error = orig - wimage;
w_sqerr = sum(sum(error.^2))

% Plot Thresholded Reconstructed Image
clf
subplot(221)
colormap(gray(256));
imagesc(gimage);
txt = ['Thresholded Reconstructed Image'];
title(txt)

% Use Haar Wavelet edge sharpening

% Pad image
pict = gimage;
pict(256,256) = 0;

imax = 8;
imin = 8;

% Call transform application
for i = imax:-1:imin
  b = 2^i;                % Get number of pixels to work on
  pict(1:b,1:b) = trans(pict(1:b,1:b));    % Work on rows of image
  pict(1:b,1:b) = trans(pict(1:b,1:b)')';  % Work on cols of image
end

% Plot input image
subplot(222)
colormap(gray(256));
imagesc(pict);
txt = ['Reconstructed Image after Haar wavelet transform'];
title(txt)

% Thin edges

subim1 = pict(1:128,129:256);
pict(1:128,129:256) = rpeak(subim1);
subim2 = pict(129:256,1:128);
pict(129:256,1:128) = rpeak(subim2')';

subplot(223)
colormap(gray(256));
imagesc(pict);
txt = ['Reconstructed Image after Haar Transform and edge filtering'];
title(txt)

% Call inverse transform application
for i = imin:imax
  b = 2^i;                % Get number of pixels to work on
  pict(1:b,1:b) = itrans(pict(1:b,1:b)')';  % Work on cols of image
  pict(1:b,1:b) = itrans(pict(1:b,1:b));    % Work on rows of image
end

subplot(224)
colormap(gray(256));
imagesc(pict);
txt = ['Reconstructed Image after Haar edge filtering'];
title(txt)


% Display squared error between original image and reconstructed image
% run through a Haar Wavelet Edge Sharpening filter

himage = pict(1:250,1:250);
error = orig - himage;
h_sqerr = sum(sum(error.^2))
